The goal of this document is to produce several interactive plots for ATAC data

# Import libraries
library(plotly)
library(heatmaply)
library(chromVAR)
library(motifmatchr)
library(BSgenome.Hsapiens.UCSC.hg19)
library(SummarizedExperiment)
library(BuenColors)
library(chromVARxx)
library(chromVARmotifs)
library(data.table)

Initialize parallel processing

library(BiocParallel)
register(MulticoreParam(4)) # adjust according to your machine

Import/Filter data

peakfile <- "../../data/fromMA/peak-ids.bed"
peaks <- getPeaks(peakfile)
counts_ma_in <- fread("../../data/fromMA/combined-counts-clean-subject-corrected-quantiles-ma.tsv")
metadata_ma <- fread("../../data/fromMA/metadata-clean-ma.tsv")
counts_mat <- round(data.matrix(data.frame(counts_ma_in)[,-60]),0)

# Make Summarized Experiment
counts <- SummarizedExperiment(
  rowRanges = peaks, 
  colData = metadata_ma,
  assays = list(counts = counts_mat)
)

Get GC content/peak; get motifs from chromVARmotifs package

counts <- addGCBias(counts, genome = BSgenome.Hsapiens.UCSC.hg19)
data("human_pwms_v2") # also human_pwms_v1 with more motifs (typically redunant)
motif_ix <- matchMotifs(human_pwms_v2, counts, genome = BSgenome.Hsapiens.UCSC.hg19)

Compute deviations; typically time consuming

dev <- computeDeviations(object = counts, annotations = motif_ix)

Find variable motifs

variabilityAll <- computeVariability(dev)
plotVariability(variabilityAll, use_plotly = TRUE)

Figure 1a; variability across motifs

Examine variable motifs in sample

mostvariable <- tail(sort(variabilityAll$variability, index.return = TRUE)$ix,30)
m <- assays(dev)[["z"]][mostvariable,]
rownames(m) <- rowData(dev)$name[mostvariable]
heatmaply(m, colors = jdb_palette("solar_extra"))

Figure 2a; motifs x samples

Plot in tSNE

#+ cache = FALSE, message=FALSE, warning=FALSE, fig.height = 10, fig.width = 10, echo = TRUE, fig.align=‘center’, fig.cap = “Figure 3a; sample tSNE in motif space”

plotDeviationsTsne(dev, deviationsTsne(dev, perplexity = 10, threshold = 2.5), sample_column = "cell.type",
                    shiny = FALSE)[[1]] %>% plotly::ggplotly()